function vmd_Fspectra(fbase, fnum)
%
%   vmd_Fspectra(fbase, fnum)
%
%   Reads file with density distribution.
%   Creates two text files with specral data.
%
%   Input:
%        fbase   base of file number (runname)
%        fnum    file number
%
%   Domain information and colormap limits are edited in the file.

%---------------------------------------------------------

% system size

%jx=48; jy=48; jz=16; jv=256; jdir=6;
jx=48; jy=1; jz=16; jv=256; jdir=2;

% physics paramters

vmax=10; Lx=450; Ly=450; 
kz=0.34; nz=1; vPhase=3.53235;


%---------------------------------------------------------

  dkz = kz/nz;
  dz = 2*pi*nz/(jz*kz);
  z  = (0:jz-1)*dz;

 
  dv = 2*vmax/jv;
  dkv = 2*pi/(jv*dv);
  v  = (0:jv-1)*dv - vmax - vPhase;
  f0 = ( 1/sqrt(2*pi) )*exp( - 0.5*(v + vPhase).^2 );


  ixc = jx/2+1;
  iyc = min(jy/2+1, jy);

%---------------------------------------------------------


  fn = num2str(fnum,'%04d');


  %-- read density distribution

  fname = [fbase, '.dfu.', fn];

  fid = fopen(fname, 'rb');
  dfu = fread(fid, jv*jx*jy*jz*jdir, 'double');
  fclose(fid);

  dfu = reshape(dfu, jv,jz,jy,jx,jdir);
  dfu = sum(dfu, 5);
  dfc = dfu(:,:,iyc,ixc);         % dfc(iv,iz)

  clear dfu
   

 
  %-- computing both spectra --

  MeanFTvSq = zeros(jv,1);

  for iz=1:jz 
    fvc = dfc(:,iz);
    fvc = fft(fvc)/(jv);
    MeanFTvSq = MeanFTvSq + abs(fvc).^2;
  end

  MeanFTzSq = zeros(1,jz);

  for iv=1:jv 
    fvc = dfc(iv,:);
    fvc = fft(fvc)/(jz);
    MeanFTzSq = MeanFTzSq + abs(fvc).^2;
  end


 %-- spectum of backgroun distribution

  ff0 = fft(f0)/(jv);
  ff0 = ff0.*conj(ff0);

 %-- write v-spectrum --

  fname = [fbase, '.DensityFTv.', fn, '.txt'];
  fid = fopen(fname, 'wt');

  fprintf(fid,'%% v-spectra in z-averaged density distribution along center line\n');
  fprintf(fid,'%% 1.kv   2.<FTvSQ>   3.F0FT^2   4.Cum(kv^4*<FTvSQ>)\n\n ')

  MySum = 0;
  for i=1:jv/2-jv/6
    k=(i-1)*dkv;
    MySum = MySum + k^4 * MeanFTvSq(i);
    fprintf(fid, "%14.6e %14.6e %14.6e %14.6e\n", k, MeanFTvSq(i), ff0(i), MySum);
  end
  fclose(fid);


  %-- write z-spectrum --

  fname = [fbase, '.DensityFTz.', fn, '.txt'];
  fid = fopen(fname, 'wt');

  fprintf(fid,'%% z-spectra in v-averaged density distribution along center line\n');
  fprintf(fid,'%% 1.kz   2.<FTzSQ>   3.Cum(kv^4*<FTkSQ>)\n\n ')

  MySum = 0;
  for i=1:jz/2-jz/6
    k=(i-1)*dkz;
    MySum = MySum + k^4 * MeanFTzSq(i);
    fprintf(fid, "%14.6e %14.6e %14.6e\n", k, MeanFTzSq(i), MySum);
  end
  fclose(fid);
    

end
